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Abstract. We study the dynamics of covariances in a chain of harmonic oscillators 
with conservative noise in contact with two stochastic Langevin heat baths. The noise 
amounts to random collisions between nearest-neighbour oscillators that exchange their 
momenta. In a recent paper, [S Lepri et al. J. Phys. A: Math. Theor. 42 (2009) 
025001], we have studied the stationary state of this system with fixed boundary 
conditions, finding analytical exact expressions for the temperature profile and the 
heat current in the thermodynamic (continuum) limit. In this paper we extend the 
analysis to the evolution of the covariance matrix and to generic boundary conditions. 
Our main purpose is to construct a hydrodynamic description of the relaxation to 
the stationary state, starting from the exact equations governing the evolution of 
the correlation matrix. We identify and adiabatically eliminate the fast variables, 
arriving at a continuity equation for the temperature profile T(y, t), complemented by 
an ordinary equation that accounts for the evolution in the bulk. Altogether, we find 
that the evolution of T(y, t) is the result of fractional diffusion. 
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1. Introduction. 

The chain of harmonic oscillators maintained out of equilibrium by means of stochastic 
heat reservoirs, is one of the few systems for which the nonequilibrium invariant 
state has been obtained rigorously pQ. However, due to the lack of any phonon 
scattering mechanism, its heat conductivity k diverges linearly with the size of the 
chain. Consequently, the Fourier's law of heat conduction Jq = —kVT, relating the 
heat flux Jq with the imposed temperature gradient VT does not hold. As a matter 
of fact, the integrability of the harmonic chain makes impossible for the system to 
support a temperature gradient. It was already recognized by Debye that the presence 
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of nonlinearities in the dynamics is a necessary condition for the occurrence of a normal 
transport, i.e. a finite heat conductivity in the thermodynamic limit. However, years 
later Fermi, Pasta and Ulam would found that nonlinear dynamics does not necessarily 
induces a statistical behaviour [2]. 

In the last decade, numerical simulations and analytic arguments have contributed 
to clarify the role of nonlinearities in the thermodynamic limit [31 HI El El 13- I n 
some studies, anharmonicities have been introduced by means of self-consistent local 
thermostats, [3, El [9]. Further attempts to derive Fourier's law in deterministic systems 
have been reported (see review papers [TUl [TTj and references therein). However, no 
rigorous derivation exists of the necessary and sufficient conditions for the validity 
of Fourier's law. Moreover, there is still a number of open questions concerning the 
steady state, such as the role of Boundary Conditions (BC in the following), while the 
convergence towards the stationary state is even less explored. 

Stochastic models are rather useful in that they can effectively reproduce the 
evolution of deterministic nonlinear systems, while allowing for analytic solutions. 
Bolster li, Rich and Visscher considered harmonic chains in which each oscillator is in 
contact with a stochastic thermal reservoir [Hj. Then, the stationary state is obtained 
assuming a self-consistent condition, namely the energy current between the local 
reservoirs and the respective oscillator is zero. Recently, it has been proved by Bonetto 
et. al. that this linear model leads to a Gaussian invariant measure and the temperature 
profiles are linear [5]. A drawback of this model is that, strictly speaking, energy 
is not conserved by the bulk dynamics (see also p2] where energy current from the 
reservoirs becomes zero in the long time limit, and [9] for their treatment in terms of 
nonequilibrium Green function formalism). Another model that can be explicitly solved 
is the Kipnis-Marchioro-Presutti (KMP) lattice model, in which stochastic collisions mix 
the energy of neighboring particles, conserving the total energy, but not the momentum 
[T3] . This model satisfies Fourier's law and a linear temperature profile is obtained. 
Energy conserving stochastic noise has been also used in lattice model systems, as 
natural generalizations of KMP and of the single exclusion process (SEP) |16j . 

In the same spirit of KMP, Basile et al. have recently studied a harmonic chain 
in the presence of random collisions among triples of nearest-neighbour oscillators |14j . 
This latter process, which conserves both energy and momentum, amounts to a diffusion 
on the energy shell. In this system, it is proven that the energy-current autocorrelation 
decays as t _1//2 and thereby heat conductivity diverges with the size of the system N, 
as k ~ N 1 / 2 [13]. More recently, Jara et al. have studied the relationship between 
anomalous heat transport and fractional diffusive processes [15j. They find that in the 
infinite system, the dynamics leading to anomalous transport is obtained from a Levy 
stable process that corresponds to a space-time scaling given by the fractional diffusive 
operator dt — V 3//2 , where V is the gradient operator. 

In a recent paper [T7], we have studied a harmonic chain with both energy- and 
momentum-conserving noise (and fixed BC). The model, a slight variant of the one 
introduced in [TJ], is amenable to analytical calculations, as the evolution equations 
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for the covariance matrix are linear (see also [IB])- Taking a suitable continuum limit, 
recalled in section [2731 we obtained a solution for the covariance matrix in the stationary 
state t — > oo. From this solution we derived exact expressions for the temperature profile 
and the heat current, finding that the heat conductivity diverges as in [13]. Remarkably, 
the temperature profile (equation (8) in [17]) is parameterfree, suggesting that it may 
represent a wider class of systems. To our knowledge, this is the first example of an 
analytic expression for the temperature profile in a system characterized by anomalous 
heat transport and confirms that a nonlinear shape persists in the thermodynamic limit. 

In this and in the companion paper [19] we extend the previous analysis to 
the nonstatonary case, obtaining a hydrodynamic equation ruling the relaxation of 
covariance matrix towards the stationary state. More precisely, here we analytically 
investigate the continuum limit, by progressively eliminating the fast variables. The 
second part contains a detailed numerical analysis of all aspects that are not (easily) 
accessible to an analytic investigation, including the case of free BC, and the estimate 
of finite-size effects. At variance with normal heat conduction, the role of BCs is very 
important in the anomalous case. For instance, in disordered chains of linear oscillators, 
the same system may even behave as a thermal superconductor or as an insulator, by 
simply switching from free to fixed BC [3]. In the present context, we find that in the 
thermodynamic limit and given bath temperatures, fixed and free BC give rise to the 
same scaling behaviour, but macroscopically different values of the heat flux. Even more 
surprising, we find that the stationary value of the heat flux varies the coupling strength 
with the heat baths only for free BC. 

The numerical analysis [19] demonstrates that these effects are caused by boundary 
layers, where the scaling behaviour of some observables changes with respect to the bulk. 
This phenomenon, which is associated with strong deviations from local equilibrium, 
hinders the development of an analytical solution in the general case. Nevertheless, for 
fixed BC, we demonstrate that the boundary layers do not affect the relevant physical 
observables, and an explicit solution can be obtained for the evolution of the temperature 
profile. Our results, obtained for large but finite chains, are consistent with those derived 
directly in the infinite-size limit [15]. This was not a priori granted in the presence of 
long-range correlations. 

This paper is organized as follows. In Section [2j we define the stochastic model and 
introduce the main notations. More precisely, in Section |2~T| we introduce the covariance 
matrix, by adopting a different definition with respect to [T7J, so as to be able to treat 
both free and fixed BC. The corresponding ordinary differential equations are derived 
in sections 12.11 and 12.21 with reference to the bulk and boundaries, respectively. In 
Section 12.31 we perform a further change of variables to simplify the treatment of the 
continuum limit, discussed in section 12.41 There, we define the smallness parameter, 
e = \/y/~N (N is the chain length) and thereby map the discrete spatial indices i, and 
j of the correlation matrices onto two continuous variables (continuum limit) x and y. 
Moreover, we introduce an Ansatz for the scaling behaviour of the different variables, as 
suggested by the numerical analysis presented in the companion paper [19]. The internal 
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consistency of the resulting equations confirms a posteriori the correctness of our initial 
choice. In section [31 we anticipate the main results to allow the reader appreciating 
them without being distracted by the technical details. The derivation of the partial 
differential equations and the elimination of the fast degrees of freedom is illustrated 
in sections HI El with reference to the bulk and boundary dynamics, respectively. The 
relaxation towards the stationary state is then discussed in Section [61 The numerical 
computation of the spectrum of the covariance evolution operator corroborates the 
analytical results. Some concluding remarks are finally presented in Section [71 



2. Stochastic model. 



We consider a homogeneous chain of N harmonic oscillators of unit mass and frequency 
oj. The first and iV-th oscillators are coupled to Langevin heat baths at different 
temperatures. The equations of motion of this system are 

Qn Pn /j\ 

p n = uJ 2 (q n +i - 2q n + g n _i) + 5„,i(£ + - Xqx) + 5 njN (£~ - Xq N ) . 

Here p n and q n are the momentum and displacement (from its equilibrium position), 
of the n-th oscillator and £ are independent Wiener processes with zero mean and 
variance 2\ksT±, where ks is the Boltzmann constant and A is the coupling constant. 
Additionally, the deterministic dynamics is perturbed by random binary collisions 
between nearest-neighbour oscillators occurring at a rate 7. These collisions are defined 
so that total momentum and energy are conserved. This type of stochastic noise is 
known in the literature as conservative noise. 

The phase-space probability density P(q,p,t) evolves according to the equation 

^ = (£ + A^)P, (2) 

The first term corresponds to the usual Liouville generator of the dynamics, acting on 
the probability density as 

with the 2N vector x = (q±, q 2 , . . . , qN,Pi,Pz, ■ ■ ■ ,Pn), and the 2N x 2N matrices a and 
d are 

-1 \ / \ 

w 2 g Ar J ' ^0 2Xk B T(r + rjs) J U 

where T is the average temperature (T + + T_) /2, 77 is the relative temperature difference 
(77 = (T + —T)/T), and 1 are the null and unit N x N matrices, 

(#t,i + 5 i>N ) , Sij = Sij (S it i - S i)N ) , (5) 



r 



<j 



and g is the negative of the discrete Laplacian 

Sij = 2Sij — 5i+ij — o~i,j+i ■ (6) 
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The stochastic collision generator is 

JV-l 

CcoiiP = lYl [ P (- • • >Pi+i>Pi> ■ ■ •) - P (- • • »Pi»Pi+i> •••)]• (7) 
i=i 

Each term in the sum expresses the probability balance for each elementary process in 
which the momenta of each pair j,j + 1 are exchanged with a rate 7. 

In this paper we consider either free or fixed boundary conditions and the choice 
adopted will be explicitly stated wherever needed. 

2.1. Covariance matrix 

The main subject of our analysis is the evolution of the covariance matrix of this 
system, namely the two-point correlation functions of the phase space variables. In 
order to develop a formalism that is both able to describe the case of fixed and free 
b.c, we consider the correlators of relative displacements and momentum, namely 
{Aqi = qi — qi-i,pi}. More precisely, we study the covariance matrix, 

where the correlation matrices y, z and v, of respective dimension (N — 1) x (N — 1), 
(N — 1) x N and NxN are defined as 

yi,j = (MAqj) , = (MiPj) , = (piPj) , (9) 

and (•) denotes the average over phase space probability distribution function P. Note 
that y and v are symmetric by definition, while z has no definite symmetry. 

In terms of the relative displacements, the boundary conditions are defined by 
imposing 

Aq 1 = and Aq N+1 = , (10) 
for free BC and 

Aq 1 = q 1 and Aq N+1 = -q N , (11) 

for fixed BC. 

Using {Aq^ instead of {q^} is necessary since the average of is not well defined 
when free BC are considered. This is at variance with [T7j, where we considered the 
correlators of the variables {qi,Pi}- However, one can recover the old formulation 
by noticing that there exists a simple mapping with the correlators Ujj = {qiqj), 
Zij = (QiPj) and Vy = (PiPj), studied in [T7j . namely 

Yi,j — Ujj — Uij+i — Ui+i,j + Ui+i,j+i j z i,j — Zjj — Zj_!j , and = Vy . (12) 

The evolution equations for c have two contributions: the dynamic contribution 
directly obtained from the equations of motion (TjQ), and the stochastic contribution that 
is evaluated upon multiplying ((Tj) by XiXj and thereby integrating over phase space. 
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We say that a given correlator in c is in the bulk of the system if the index of the 
momentum variable is in [2, N — 1] and the index of Aq is in [3, N — 1] for free BC and 
in [2, N] for fixed BC. The evolution equations c are in the bulk 

y*j Z i>* 1 ' ^ i , j — 1 j 

Zij = v m - Vi-i,j + cj 2 (yij+i - y-jj) + 7 ( z i,j+i + z ij-i ~ , (13) 
Vjj = ^ 2 (zj+i^ - Zj ti + Zj + ij - Zjj) + 7Wij , 

where the N x N collision matrix VF, corresponding to the contribution from the 
stochastic noise, depends on the distance between the evaluated indexes i and j and can 
be written in a compact form as 

Wij = i(A\v ; . :j + 5 jt iVij-i) + • i('\.:v, i. ; + 5j )JV Vij+i)] 

+^i,j(^i,N v i+l,j+l + <5j,lVi_ij_i) — (2(5jj_i + Sij+i — 5ij) 

- 0~i,N - - 5j,N + ^,1^3,1 + Si.NSj^Vij , (14) 

where <5jj is the Kronecker delta function and Si j = 1 — VFjj also holds for the 
boundary terms, and since it deals with momentum variables only, it is independent of 
the specific BC. 

2.2. Boundary conditions 

As a consequence of the physical boundary conditions imposed on the oscillators of 
the chain edges, the border terms of the covariance matrices, follow a dynamics that is 
different from ( TIB"]) . The BC affect the phase space variables Aq iy according to ( fTUj) and 
(jTTj) . and Pi, through the boundary terms in (fl4|) . Furthermore, the coupling between 
the oscillators at the edges and the heat baths affects the evolution of their momentum 
(CQ). In this section we content ourselves writing down the equations of motion of the 
border covariances corresponding to the first column and last row of each matrix. The 
equations for the matrix elements in the last column and first row can be obtained 
analogously. The latter are not a mirror image of the former. Nevertheless, one can 
verify that to leading order in the continuum limit, they lead to the same behaviour. 

For fixed BC the equations of motion for the border matrix elements of the first 
matrix column (index 1 for both phase space variables) are 

= <5i,Ar+i z i,i - <5i,iz M _i + z i; i , (15) 

z i,l = 5i,JV+lVi,l - ^iVi-1,1 + W 2 (yi,2 - yt,l) + 7( z i,2 - Zt,l) - , (16) 

v i)X = u 2 (z i+lil - z M + z 2 ,; - z M ) - 1(1 + 6 ifl + 5 i:N )v i:1 + S itl 2lk B T + + jW iA , (17) 

where 1 < % < N + 1 in IflSll . ffTHjl and 1 < % < N in (TIT]). For fixed BC the last matrix 
row different from zero corresponds to index N + 1 for Aq variables and to index N for 
p variables. The corresponding equations of motion are 

YN+lJ = ~ z j,N + 0~j,N+l z N+l,j — <5j,lZjV+lj'-l , (18) 

zjv+i,j = — vjvj + u; 2 (yAr + ij + i — yjv+ij) + l(o~j,N z N+i,j+i + $j,i z N+i,j-i 
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-(&j,N + 8j,l) z N+l,j) , (19) 
VJVJ = ^ 2 ( z N+l,j — z N,j + z j+l,N — z j,N) — 1(1 + <5j,l + 0~j,N)vN,j 

+6 jtN 2lk B T_ + t Wjv,j , (20) 

where 1 < j < N + 1 in flU} and 1 < j < N in (JTSJ), fl20|). 

For free BC the first matrix column (different from zero) is index 2 for Aq variables 
and 1 for p variables. Consequently, the equations of motion are 

ji,2 = z 2,i - Z 2) i-1 + z i.2 - Zj ; i , (21) 

Zi,i = Vi,i - Vi_i,i + ^ 2 yi,2 + 7( z i,2 - Zi,i) - izi,i , (22) 
Vj,! = ^ 2 (<5 iiA rz i+li i - 5uZi,i + z 2) j) - 1(1 + + 4tv) v o + 8 i , 1 2lk B T + + , (23) 

where 2 < i < N in (|2ip . (1221) and 1 < z < A" in (1231) . For the border matrix elements 
of the last row (index A" for both Aq and p variables), the equations of motion are 

yjV,j = z j,N — z j,N-l + z N,j — ZjVj-1 , (24) 
+"f(8j,N z N,j+l + ^j,lZjV,j-l - (^',1 + $j,N) z N,j) , (25) 

vjv,j = ^ 2 ( - zjvj - z j,N + 8 jtN z j+ltN ) - 1(1 + S jtl + 5 j:N )v N:j 

+6 j>N 2lk B T_ + ^Wnj , (26) 

where 2 < j < iV in (ED and 1 < j < N in ([25]), (1261). 

The equations presented in this section (together with the equations for the border 
elements of the first row and of the last column of the matrices) constitute the dynamic 
boundary conditions of (1T31) . 

2.3. Change of variables 

We shall see in the following sections that certain combinations of covariances appear 
naturally in the evolution equations. It is therefore convenient to introduce some of 
these combinations as new variables. An important example of such combinations is 

^ij = V M - ^Vtj , (27) 

which has a precise physical meaning: its diagonal elements ipi^ correspond to the local 
balance between kinetic and potential energies. If our system satisfies the virial theorem, 
then we should find that ip^i = for all i. As we will discuss in the following section 
15.11 this is not always the case. In the rest of the paper we will substitute y in favour 
of 

In jTTj we found that the stationary state solution of (1T31) obtained by taking all time 
derivatives to zero, implies that the covariance matrix Z is antisymmetric. From (1T2|) it 
is clear that the covariance z has no definite symmetry, not even in the stationary state. 
Therefore, it is pertinent, as we do in the following, to decompose z in its symmetric 
and antisymmetric components z 1 * 1 

z^z+. + zr., with z± (28) 
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= u 2 [ - 4z+ + z+ +1 + z+ 1;j + zf_ l;j + zf^ - zr j+1 + z r +1) . - z r 1J + z r ? _ 1 ] 

+7^ , (29a) 

= 7 [ Z i,j+1 + Z i+l,j + Z i-lj + Z i,j-l + Z i'+1 ~~ Z i+lj ~~ Z i-lj Z ij-1 ~~ ^Zjj] 

- ^ij+i + Vi,j+i - Vj+ij - Vi_ij + Vij-i , (296) 

2 z ij = 7 [ z ij+l + Z i+lj + Z t-l,j + Z i,i— 1 + Z i,j+1 ~ Z i+l,j ~ Z i-l,j + Z ij-1 ~~ ^ Z i'] 

+2^ij - - V'ij+i + v i,i+i + v i+i,i - Vj-i,j - Vij-i , (29c) 

= ^ [-2z+- + + z ti,j - z w + z T+i,j] + -rWij , (29d) 
Other useful combinations of covariances will be considered when needed. 



2.4- Perturbative expansion and continuum limit 

Our first goal is to transform the set of difference equations (129 aft - (129 dp into a set of 
partial differential equations, by taking a continuum limit that is appropriate to our 
problem. We do this by following a perturbative-like analysis and choose e = l/y/N 
as perturbation parameter. In order to proceed forward we need first to attribute the 
right order (in powers of e) to the covariance matrix elements. In order to keep the 
presentation as simple as possible, we proceed on the basis of our knowledge of the 
stationary solution [17] and of the numerical solutions discussed in [19] . 

In the stationary state, y id and Vjj are 0(e), with the exception of the diagonal 
terms that are 0(1), being proportional to the mean potential and kinetic energy, 
respectively. On the other hand the combination in equation ( J27j) is 0(e 2 ), indicating 
that the system is locally at equilibrium. Moreover, it can be shown that the fields z ± 
are x-derivatives of Z, namely z + oc Z x and z~ oc Z xx . In jTTJ we found that off-diagonal 
terms of Z are 0(1), while along the diagonal Z = due, to its antisymmetry. Since 
each differentiation w.r.t. x increases the order by e (see ( 1321) below), z + must be 0(e) 
and z~ of 0(e 2 ). Altogether, 

ipi,j = e 2 5 itj tt i + e 2 {l-5 iJ )y itj , 

z tj = eSijS? +e(l- 6 itj ) Z+ , 

z 7, = e 2 (1 - 5 hj ) Z^ , 

v i,j = fiijTi + e (1 — 5ij) Vij . 
These equations fix also the notation that we will use in the rest of the paper. Note 
that in this notation, all the correlators appearing at the r.h.s. of the definitions above, 
namely Q, S + , Z + , Z~, T and V, are 0(1). We refer the reader to [19], where we give 
numerical support for the validity of (130|) . Furthermore, we have found in [19] that for 
fixed BC, the off-diagonal elements of i/j scale as e 3 for A = u, and as in ( 1301) for A^w. 
In what follows we will assume (1301) and discuss the modifications for the resonant case 
A = uj where appropriate. 

The last step before proceeding to the derivation of the partial differential equations 
consists in transforming the discrete indexes of the correlators into two suitable 
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continuous variables. We do this by introducing the variable y for the diagonal direction 
and variable x for the transversal direction as (see e.g., figure 1 in [TT]). 

c -\ (i+j)e 2 -l 
x = (i-j)e; y = J^J-—^- ( 31 ) 

The nonlinear definition of y is chosen so that its domain [—1, 1] is independent of x. 

Nevertheless, in the limit N — > oo the effects of the nonlinearities are localized at the 

boundaries of the domain and thus, do not complicate the study of the bulk dynamics. 

Differential changes in these variables can be written as 

, (i + j)e 2 - 1 + se 2 

x = x + fe ; y = - — — , 32 

1 — \i — j\e 2 — fe 2 

where the integer shift functions f(Ai,Aj) = Ai — Aj and s(Ai,Aj) = Ai + Aj, 
(/, s : Z 2 I— * Z), account for the displacements on the discrete variables along the x and 
y directions respectively. 

Using this, the continuum limit of any covariance matrix element irij+Aij+Aj can 
be written up to 0(e 2 ) as 

m i+A i,j+Aj = m(x + fe,y + e 2 (fy + s)) , for i > j , (33) 

for the continuous correlator function m. Here and in what follows, we keep the 
same notation for the continuous functions derived from the matrix variables (e.g. 
v i,j{t) — > v ( x iy,t) etc.). If we restrict to the matrix lower triangle, as we do in 
( 1331) and in the following, then in the limit N — ► oo the variables (x, y) belong to the 
domain 

V = {(x,y)\xe [0,oo); yE [-1,1]} (34) 

Note, for instance, that x = const corresponds to moving along the diagonal direction, 
{x = corresponding to the main diagonal). For more details we refer the reader to 
section 4.1 of [TT] . 

3. Main results 

In the next section we show that, after eliminating the fast degrees of motion, the 
bulk dynamics reduces to the following two equations (the subscript denoting partial 
derivation) 

= e 2 (jZ+ x + 2V y ) , (35) 

V = e 2 (jV xx + 2u 2 Z+) , (36) 
together with the conditions on the boundary of the domain 

5V = {x = 0} U {y = ±1} . 

On x = 0, we find that 



V x (0,y,t) = 0, (37) 
Z+(0,y,t) = -^T y , (38) 
f(y,t) = 2e 3 uj 2 Z+(0,y,t) . (39) 
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This last equation is accompanied by its BC, namely T(=pl,t) = T±. The peculiar 
mathematical structure of the problem should be underlined: the boundary evolution of 
Z is determined dynamically by the differential equations ( 1381) and ( 1391) which, in turn 
determine the temperature field T(y,t). 

On y — ±1 the BC are particularly difficult because close to these boundaries, the 
expansions fl30l) are no longer valid. The numerical solution in [19], shows that a kind of 
"boundary layers" (BL) exist, namely that at least some of the fields scale differently in 
the regions of T> lying within a distance 0(e) from y = ±1. In particular, in this region, 
ijj is of order 0(e). Recalling that ip is the difference between the kinetic and potential 
energy, which away from the diagonal are 0(e), this implies that in the BL the relative 
difference between potential and kinetic terms is of 0(1), which means that the system 
is not in local thermal equilibrium. 

The existence of a BL hinders us from finding an explicit exact solution in the 
general case. This technical difficulty is irrelevant in two cases: a) for fixed BC and b) 
for free BC at the "resonant" value of the bath coupling constant, A = oj. In case (a), 
the boundary layer does not affect the relevant fields and, as we show in section [51 it is 
sufficient to impose 

Z + (x,±l,t) = , V(x,±l,t) = . (40) 

This allows us to find an explicit time-dependent solution. From a physical point of 
view, case (b) corresponds to the only situation in which the coupling can be perfectly 
tuned to avoid an impedance mismatch on the boundaries. Unfortunately, we were not 
able to find an explicit solution in this case. 

Inspection of the equations of motion reveals that the temperature field T is the 
slowest variable (it evolves on a time scale 0(e~ 3 )). Since Z + and V relax on a time 
scale 0(e~ 2 ) they can be adiabatically eliminated by setting their time derivative equal 
to zero in equations fl35|) and fl36|) . By further eliminating the field V, we obtain the 
fourth-order equation 

j 2 Z+ xxx - 4uj 2 Z+ = , (41) 

that is formally equal to the equation solved in [IT], the main difference being that here 
Z + is time dependent. The dynamical equation is obtained by first determining the 
stationary solution of fj4"Tj) with the appropriate BC and then using fl38l) to express Z + 
as a function of T and replacing the result into ( |39l) . This is accomplished by considering 
the Fourier expansion 

oo 

T(y,t) =T s (y) + Y,Tn(t) sin 

n=l 

with T s being the stationary solution of T, as given by formulas (18) and (19) of [17] . 
In section [6] we show that the coefficients T n (t) obey the linear equation (1741) . The 
associated eigenvalues, that must be computed numerically as the problem is not exactly 
diagonal, uniquely determine the relaxation to the steady state for any assigned initial 
condition. They are found to be proportional to — (k/N) 3 ^ 2 (for k being a positive 



nix , 



(42) 
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integer labeling the eigenvalues). This is reminiscent of the spectrum of the eigenvalues 
of the fractional Laplacian V 3 / 2 , thus suggesting that the evolution of the temperature 
field is ruled by some underlying fractional diffusion equation on hydrodynamic scales 
[151 EDI. 



4. Dynamical equations 

We first focus on the dynamics in the "bulk" of system, namely the interior of T>. In 



Appendix A we derive the set of coupled partial differential equations ( 11.3al) - (11.3aT) . for 
the covariances ip, z~, z + and v. Furthermore, using (133]) the equations with their 
explicit order in e are rewritten as 

V = 2e {2u 2 Z x + ou 2 Z+ x + jV xx ) + te 2 uu 2 yZy , (43 a) 
ST = eV x + e 2 (>yZ- + y%) , (436) 

Z + = e 2 tfZ+ + 2V y ) - e 3 (U> xx + tf v ) , (43c) 

V = e 2 (2u 2 Z- + uu 2 Z+ x + 2^V XX + 2uu 2 Z+) + 2e 3 u 2 yZ~ . (43d) 

The four variables can be split into a pair of fast (\1/ and Z~) and slow {Z + and V) 
ones, which evolve on time scales of order e~ x and e~ 2 , respectively. Upon substituting 
the x-derivative of (143 a\i into the time derivative of (143 61) we find 

Z- - As 2 uj 2 Z xx - 2e 2 (uj 2 Z+ xx + jV xxx ) = 

5 2 7 Z" + 8e 3 u 2 yZ xy + 2e 3 y {u 2 Z+ y + V xxxy ) , (44) 

where we have retained only terms up to order e 3 . The terms on the r.h.s. do not affect 
the final solution but must be taken into account to justify the adiabatic elimination. 
In fact, scaling the time by e, we see that they are o(e 3 ) and could, in principle, be 
neglected. However, if we do so, we are left with a non-dissipative wave equation (with 
source terms) for Z~ , that cannot account for the convergence towards the steady 
state. Therefore, to study the evolution of the fast variables, we are obliged to include 
the higher order terms (losses are actually provided by the Z~ x term, the other being 
perturbations of the source term). After this remark, we are authorized to adiabatically 
eliminate Z~ and From (143 al 143 61) we obtain to leading order 

Z- = - \zt x + ^V xx , (45a) 

V x = 0. (456) 

Note that \& is constant moving away from the diagonal. We now turn our attention to 
the slow variables and substitute the above two equations into (143 c\ 143 dj) . To leading 
order, we finally obtain (1331) and (133]) . 

At this point, it is useful to illustrate some features of (1331) and (133]) . The symmetry 
of these equations suggests to introduce the new variables = ujZ + ± V, that allow 
decoupling the system of equations into 

QW = e 2 ( 7 QW±20 . (46) 
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The first term on the r.h.s. describes a transversal diffusion process characterized by 
the diffusion constant e 2 7; the second term accounts for a longitudinal right/left sound- 
wave propagation, depending on whether Q~ / Q + is considered. Leaving aside for the 
moment the issue of BC, by absorbing the order e 2 in the time variable, we look for 
solutions of the form 

Q + {x,y,t) = P(x,t) cos K y y + Q(x,t) sin K y y , (47) 

we consider only Q + as the equation for Q~ leads to the same dispersion relation). By 
substituting this in ( 146|) and separating the independent terms, we obtain 

P = jP xx + 2uK y Q , (48) 

Q = iQxx - 2uK y P . (49) 

By then assuming the following form of the solution, 

P(x, t) = a(t) cos K x x + b{t) sin K x x , (50) 
Q(x, t) = c(t) cos K x x + d(t) sin K x x , (51) 

we find that a(t) and c(t) satisfy a system of two ordinary differential equations, 

a = - ^K 2 a - 2uK y c , (52) 

c = - -fK 2 c + 2uK y a , (53) 

(b and d do not add any information as they satisfy the same set of equations). Looking 
for solutions of the form a(t) = aexp(/xt), c(t) = cexp(/xi), the resulting eigenvalue 
equation yields two degenerate branches for the dispersion relations, 

H = -ryKl ± 2uoiK y , (54) 

where % denotes the imaginary unit. 

Taking \i = 0, we recover the eigenvalue of the stationary state [T7]. Most 
important, the real part of \i is negative, thus ensuring the stability of the stationary 
state. By reintroducing the e 2 factor in the time units, we can thus conclude that 
modes characterized by a K x of 0(1) relax on a time scale of order e~ 2 . As discussed 
in [19], these are the modes that mostly contribute to the relevant nonzero off-diagonal 
correlations. However, K x can, by construction, be as small as e~ x . As a result the 
slowest relaxation times that one can observe are of the order of iV 2 . This is indeed 
confirmed by the numerical calculation of the spectrum of the evolution operator |19j . 



5. Boundary conditions 

A complete solution of the dynamical problem requires solving the set of partial 
differential equations (143aj) - (143d!p on the whole domain T> (1341) . including its boundary 
at all times. In a general context, the difficulty of obtaining a full solution depends on 
the constraints along 8T>. If they amount to algebraic conditions or if the boundary 
dynamics is faster than the bulk dynamics, then one deals with standard type of static 
BC. In the opposite case, it is the bulk dynamics that can be adiabatically eliminated 
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and the relevant (long-term) evolution would be controlled by what happens along the 
boundaries. 

In the first section we study the dynamics of the covariance elements close to the 
diagonal and derive a differential equation describing the evolution of the temperature 
profile. Later, in section 15.21 we study how the physical BC determine the dynamics of 
the correlators along the boundaries y — ±1. 



5.1. Boundary conditions along x= 

The BC along the diagonal are not connected with the physical BC of the chain, 



but rather with the symmetry of the various matrices. In Appendix B we derive the 
differential equations describing the evolution of the covariances at x — 0, f l2.3fflH2.3cj) 
and fl2.5aH2.5dj) . As the diagonal terms Q and S + appear only as differences with their 
respective off-diagonal counterparts, it is convenient to introduce 5^ = ^ — Q and 
s5Z + (y,t) = Z + (0,y,t) — S + (y,t), where we benefit of the numerical investigations to 
anticipate that the latter difference is of higher order than the single addendajjjl 

By introducing 5^ and 5Z + and using (l30j) . we obtain the following set of partial 
differential equations describing the evolution at x = 0: 

5V = -2u 2 (35Z + + Z~ + 2Z+) + 2 7 V X + e( - 4uo 2 yZ+ + 7 /C v ) , (55a) 
^ = 2u 2 (Z- - 5Z + ) + 27V,. + e(2u 2 (Z+ x + 2Z~) + 7 /C v ) , (556) 

Z~ = ^ - iZ~ + T y + e(tt x - V y ) , (55c) 
5Z + = ^ - 1 (35Z + + 2Z+) +T y + e(V x - 2jyZ+ - V y ) , (55 d) 

Z + = e(±89 - !5Z + + T y ) + e 2 (V y + >yZ+) , (55c) 

V = e(u 2 (Z~ - 5Z + ) + 2 7 V X .) + e 2 (uj 2 (Zt x + 2Z~ + 2Z+) + 7 /C v ) , (55/) 
f = 2e 2 u 2 (5Z + + Z- + Z+) + s 3 u 2 (Z+ x + 2Z X + 2(y + \)Z+) , (55a) 

where we have considered the first two leading contributions to the evolution of each 
variable and, for the sake of compactness, we have introduced 

K.v = 3V ra + 2yV y . (56) 

Like in the bulk, the variables evolve over manifestly different time scales, and the 
temperature field T is the slowest one. Therefore we proceed by adiabatically eliminating 
the other variables, starting from the fastest ones. By setting the first three time 
derivatives equal to zero and considering only the leading terms, equations (155 aj) . ( 15561 ) 
and (155 cl) lead to 

^ = -(7^ + + ^V, + 27;) , (57a) 

§ As a matter of fact, assigning to 8Z + the same scaling as its addenda leads to unphysical super fast 
evolution. 
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sz+ = 5 B v * - ^ + ) ■ < 57c > 

Moreover, the stationary solution of (155 cfl) yields 

7^ + ^V, + T y = 0. (58) 

By now inserting 1 15 7 d) and ( 15 7 eft into ( 155 eft , and setting the time derivative 2 + equal 
to zero, we obtain (1371) . This equation has an obvious meaning: V is symmetric by 
definition across the diagonal, so that we naturally expect V to have be maximal for 
x = 0. 

Furthermore, by using (1371) in (158|) . we obtain the second relevant constraint (138|) . 
From ( 15761) and ( 15 7 eft , one can easily find that 2~ = 5Z + . By then going back to (I57a[) 
( 157&P and (I57c|) . and using the constraints (1371) and (|38|) ). we obtain 

6Z + = Z- = j-T y , <J* = -T y . (59) 

Altogether, once the five conditions contained in ( 1371) . ( 1551) and ( 1591 are satisfied, it turns 
out that the derivatives of all seven boundary variables (l55all55gD are equal to zero to 
leading order. In particular we see that the variable ^ remains undetermined, but this 
is not a problem, as it is of higher order in e and does not contribute to the leading-order 
evolution of the physically relevant variables. Additionally, with the exception of V, all 
the variables are expressed as a function of the temperature profile T, whose evolution 
must be determined if we want to find a closed solution. 

It turns out that the leading contribution of 0(e 2 ) to T is zero. Therefore, it is 
necessary to go one order further in the perturbative analysis. This can be easily done 
by noting that the leading contribution to T is equal to that of Cl — ^ — 5^ (see 12. 3d) . 
Subtracting ( 155 6h from ( 155 ah and setting the time derivatives equal to zero we find that, 
up to 0(e 2 ), 

2u 2 (5Z + + Z- + Z+) = -slu 2 (Z+ + 2Z- + 2yZ+) . (60) 

By inserting this into ( |55gP and retaining terms up to 0(e 3 ), we obtain equation (1591) . By 
recalling that Z+ is proportional to the divergence of the heat flux along the diagonal, we 
recognize that (1591 is nothing but the continuity equation for the energy and could have 
been derived simply on the basis of physical arguments. However, the relaxation of the 
temperature profile towards the stationary state, occurs on a time scale that is 0(e~ 3 ), 
i.e., for t ~ N 3 / 2 . As a consequence, we can conclude that the bulk dynamics is faster 
than that occurring along the diagonal and can, thereby, be adiabatically eliminated as 
anticipated in section [5J 

Finally, in order to complete the treatment, we must complement (1591 with its 
physical BC, as it is a (one-dimensional) partial differential equation. Without the need 
of a formal treatment, it is easily understood that these BC are simply T(±l,£) = T T . 
As a matter of fact, the relaxation on the boundaries occurs on a finite time scale 
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(~ 1/A), i.e. it is basically instantaneous with respect to the above mentioned time 
scales. 

5.2. Boundary conditions along y = ±1 

In this section we analyse the conditions that the physical boundary conditions, either 
fixed or free, impose on the covariances. At the chain edges, where the system is directly 
coupled to the heat baths stochastic evolution, the dynamics is different from that in the 
bulk: on the one hand the deterministic restoring force is not counterbalanced by the 
boundary and on the other hand, the stochastic collision for the edge oscillators is also a 
"one-sided" process. Consequently, we introduce new auxiliary variables to distinguish 
the boundary dynamics from its bulk counterpart. More precisely, we define </>j.i = i/j^i, 
Cj = z ti an d = Vi,i. 



In Appendix C we derive the partial differential equations describing the dynamics 
of the covariance at y = —1 for fixed and free BC. However, in this case we are not 
entitled to use (|30|) in order to assign the correct order in e of the covariance variables. 
As we have discussed in Section El there exist a boundary layer, namely a region around 
y — ±1 of size e~ x , where the scaling of the covariance matrix on e differs from (130]) . 
This BL has been further studied numerically in [19]. It is important to note that if 
we insist using (130]) then mathematical consistency requires that the order of e.g., ip is 
0(e) and not 0(e 2 ), which is also what we numerically observe in [19J. However, (130]) 
cannot differentiate the scalings in the BL from those in the bulk. Consequently, in this 
Section we do not use the expansions (130]) and limit ourselves to extract some physical 
information from the leading contribution determined only by the differential structure 
of the equations. First we focus on fixed BC. 

For fixed BC we end up with four equations for the bulk variables and four equations 
for the boundary variables, (13. 3 aft - (13. 3 /?1) . Since all equations evolve on time scales of 
0(1), they can be adiabatically eliminated from the bulk dynamics, which is at least 
0(e). By taking all time derivatives to zero and solving the resulting set of relations, 
we find that in the stationary state all boundary variables coincide with the respective 
bulk counterparts, 

( + = z + , (~ = z~ , v = v and = ip , (61) 

and that 

v(x, -l,t) = 0, (62) 
z + (x, -l,t) = -z-(x, -l,t) . (63) 

It is interesting to point out that these relations are independent of the parameters of the 
system and that (]6T]) imply that the bulk variables are continuous at the boundaries. 
Recalling that in ( !63l z is of higher order than its symmetric counterpart z + , we 
do obtain ( 140]) . It is interesting to note that the indetermination of ip(x, — 1) and 
z~(x, — 1), on which the effects of the BL are mostly observed, does not affect the 
leading order dynamical solution of the physically relevant fields. Therefore, for fixed 
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BC, the existence of a BL does not impede us from using the boundary relations (j40j) 
and determining the evolution of v and z + . 

We now turn our attention to free BC. As seen in Appendix C.2 , in this case, only 
two auxiliary variables along the boundaries are necessary, = Zj ; i corresponding to 
the non symmetrizable term, and z/^i = v^i. From fl3.5ap - p.5dl) we find that when the 
three conditions 

v — v , C = z + and z~ = (64) 

are satisfied, the four time derivatives are equal to zero (to leading order), thus satisfying 
the stationary state solution. By using these relations in ( l3.5d) -( 3~5^ ), we obtain the 
relevant mathematical conditions, 

ip(x, -1, t) = v(x, -1, t) - lz + (x, -l,t) , (65) 

u 2 z + (x, -l,t) = lv{x, -l,t) . (66) 

There are two main differences between ( 1621) . ( 163]) and the equations above: first, ( |65i) . 
( 166]) depend on the variable ip, which, as we have discussed, it is the variable whose 
behaviour is most affected by the BL. Second and more important, the free BC relations 
now depend on the parameters A and uj. 

By combining (1651 and (166]) we find that 

(uj 2 - l 2 )z + (x, -1, t) = ty(x, -1, t) . (67) 

If uj — Z, then tj}(x, —1) = (or at least 0(e 2 )), consistently with our expectations 
from the bulk dynamics. In this resonant case the boundary condition reduces to 
ujZ + (x, — 1, t) = V(x, —l,t). Though simple, we have not found a way to derive an 
explicit solution of the bulk equation which satisfies this constraint. 

In the non resonant regime u 2 ^ I 2 , (l6Tj) implies that z + and ip are of the same order 
along the boundary. The numerical studies presented in [19] confirm this prediction, but 
show also that this is because ip is of 0(e). Such an observation is seemingly inconsistent 
with the bulk analysis which predicts x/i to be of higher order. As said, the only way 
to solve the paradox is by invoking the presence of a BL connecting the two different 
scaling regimes. 

Summarizing this section, the BC for y = ±1 reveal a crucial difference between 
fixed BC and free BC. In the former case they are independent of parameters of the 
system, which in turns imply that the asymptotic profile, as well as the leading term of 
the heat flux are universal. In the latter case both quantities depend on the coupling 
strength with the heat baths I. In this very atypical situation the contact with the baths 
may lead to measurable macroscopic effects. Note in particular, that the heat flux for a 
system with anomalous thermal conductivity, like the one that concerns us, in general 
may depend on the type of boundary conditions, even in the infinite volume limit. 



6. Dynamics of the temperature field 



The evolution of the temperature profile is determined by (1591 . subjected to conditions 
(1371) and (l38l) . Since the temperature field evolves on a slower time scale, the bulk 
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dynamics can be adiabatically eliminated. Therefore, it suffices to solve (14 ip for Z + 
and plug its solution into fl39|) . As discussed above, an explicit calculation is feasible 
only for fixed BC, where the first of conditions ( f40l) implies that we can, following j!7j . 
expand the time-dependent solution of Z + as 

z + {x, y,t)=Yl B ^ f ) sin [y- (v + !)] • ( 68 ) 

n 

The Fourier coefficients satisfy the ordinary differential equation 



2 



whose explicit solution yields 

B n (x, t) = A n (t) exp(— a n x) sin(a n x) + A' n (t) exp(— a n x) cos(a n x) , (70) 

where a n = \J mxu/2^f and we have discarded the components which diverge for x — > oo. 
By differentiating the equilibrium solution of fl35|) with respect to x, we realize that the 
condition ( 137|) is equivalent to Z xxx (0,y,t) = 0, which in turns implies that A n = —A' n . 

If one is interested only in the stationary solution, equation (1391) implies that 
Z + (0,y) = constant, namely that the heat flux is constant along the chain. This 
condition transforms itself into distinct equations for the coefficients {A n }, which can 
therefore be determined (apart from a multiplicative factor). Afterwards, with the help 
of ( 1381) we can determine T s (y). The unknown multiplicative and additive factors are 
eventually removed by imposing T s (=pl) = T±. We do not report these calculations as 
they would closely follow what already reported in [T7] . 

Here we wish to solve the dynamical problem, particularly for the temperature field 
T(y, t). Let us consider its Fourier expansion ( T42l) where we have only included the terms 
that are appropriate for fixed BC. To write down closed equations for the coefficients T n , 
we must face the problem that the two sides of equation (139!) are expanded in a different 
set of functions, namely sines and cosines, respectively and the problem is therefore not 
diagonal. By using vector notations with an obvious meaning of the symbols, we obtain 
from (EHD and (E9D, 



where 



A = — DRT (71) 
f = 2e 3 u 2 RA (72) 



2k 2 /(k 2 - n 2 ) for k + n odd 
otherwise 



and D nym = 5 nm /a n . We can thus write a closed equation for T, 



£ 3 u; 2 



T = RDRT (74) 

7 

A numerical evaluation reveals that RDR is almost diagonal. In fact, the 

eigenvectors are very close to Fourier sine-modes [Hj. In figure Ofl we report the 
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eigenvalues in ascending order, versus the index I that is equal (within a proportionality 
factor) to the corresponding wave number. The data align almost perfectly along a 
straight line (in log- log scales) that corresponds to a scaling with a power 3/2. The 
deviations observed at small wave numbers are not due to the truncation of the operator 
RDR; they express the fact that RDR is intrinsically defined on a finite domain. In 
[T9] , we show that the numerical solution of the entire dynamical operator (without any 
approximation) confirms our analytical predictions. 

Altogether, the space-time scaling of ( 1711) indicates that the evolution of the 
temperature field T(y,t) is, on the considered time-scales, ruled by a diffusion equation 
with a fractional Laplacian V 3,/2 . Recently, this has been shown to be the case for a 

• i ii j |i r| t ii • ii n ' i Tvri* ' i / • *ii i i i* < i 

simi 

effec 




Figure 1. Spectrum {A^} of the linear equation ([71)) . The eigenvalues are expressed 
in e 3 uj 2 /j time units. The straight line corresponds to a power law with a rate 3/2. 



7. Discussion and conclusions 

We have presented a detailed description of the relaxation towards the nonequilibrium 
steady state in a model of harmonic oscillators with conservative noise. To our 
knowledge, this is an almost unique instance where relaxation phenomena can be studied 
in great detail in a realistic setup. By implementing the continuum-limit ideas previously 
introduced in [T7] we have obtained a set of partial differential equations describing the 
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evolution of the covariance matrix. In the bulk, the velocity-velocity and the symmetric 
component of the velocity-position correlations are the relevant (slow) variables: they 
appear to evolve on a time scale of order 1/e 2 = N. This means that in the bulk, 
relaxation phenomena are mostly controlled by the propagation of sound waves. 

Along the boundaries, the evolution of the relevant two-point correlators can be 
explicitly determined to lowest order in e. Again, these correlators evolve on different 
time scales, the temperature field being the slowest one (its dynamical equation evolves 
on a time scales t ~ N 3 ^ 2 ). By adiabatically eliminating the fast variables, we 
find that T(y, t) satisfies an energy continuity equation, where the expression for the 
current can be obtained from the stationary solution of the bulk equation. Altogether, 
the temperature field T(y,t) appears to satisfy a diffusion equation with a fractional 
Laplacian V 3 ^ 2 . However, the relationship with fractional Brownian motion should be 
further explored. In particular, it is not yet clear to what extent the temperature profile 
can be obtained from the solution of the fractional equation in a finite domain. 

The case of free BC remains open in view of the difficulties arising along the 
boundaries, where the mathematical conditions depend explicitly on system parameters 
such as uj and A. This is not only the indication of a lack of "universality" but implies 
also that some variables (namely ip) must scale differently in the bulk and along the 
boundaries. As a result, boundary layers are expected to arise (and this is confirmed 
by the numerical analysis carried out in the companion paper [19]) which would require 
a separate analysis. This is one of the open problems that will be worth investigating 
in the future, especially in the perspective that a similar scenario might hold in generic 
nonlinear deterministic systems. Only in the resonant case u 2 = A 2 , boundary layers 
do not exist. However, even in this limit, the mathematical conditions holding on the 
boundaries are sufficiently complicate to prevent the derivation of explicit expressions 
(at least, to the best of our knowledge). 

All of our analysis has been restricted to two-point correlators. The main reason 
is that the dynamical equations are exactly closed onto themselves, so that there is no 
need to invoke higher-order correlators. Moreover, this analysis allows determining exact 
expressions for the most relevant variables such as the heat flux and the temperature 
profile. However, one should not forget that the scaling behaviour of heat conductivity 
in this stochastic model (k ~ A 1 / 2 ) differs from that of generic nonlinear systems, 
where k ~ A 1 / 3 . Is this an indication that a faithful description of such systems needs 
including higher-order correlators? More modestly, it would be already interesting to 
check to what extent a Gaussian approximation of the invariant measure based on the 
knowledge of two-point correlators can accurately describe other variables such as, e.g., 
energy fluctuations. 

Appendix A. Derivation of the differential equations: bulk dynamics 

In this appendix we derive the set of partial differential equations describing the 
dynamics of the covariance matrix elements in the bulk of the system, namely for 
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covariances m^- for which 

(i) the index of the momentum variable is in [2, N — 1], 

(ii) the index of Aq is in [3, N — 1] for free BC and in [2, N] for fixed BC, 

(iii) and \i — j\ > 1. 

In this situation the stochastic collision matrix is simply 



In the continuum limit, discrete shifts of the indexes i and j yield infinitesimal 
changes of the field variables x and y and by using (|33|) . the set of difference equations 
(I29ap -( l29dl) lead to a set of continuous equation given explicitly by 



ip (x, y) = u 2 [- 4z + (x, y) + z + (x - s, y - e 2 (y - 1)) + z + (x + e, y + e 2 {y + 1)) 
+z+ (x - e, y - e 2 (y + 1)) + z + (x + e, y + e 2 {y - 1)) - z" (x - e, y - e 2 {y - 1)) 
+z~ (x + e,y + e 2 (y + 1)) - z~ (x - s, y - e 2 (y + 1)) + z" (x + s, y + e 2 (y - 1)) ] 
+7 [ v (x + e, y + e 2 {y + 1)) + v (x - e, y - e 2 {y - 1)) + v (x - e, y - e 2 {y + 1)) 



2zT (x, y) = 7 [ - 4z~ (x, ?/) + z~ (x - e, y - e 2 {y - 1)) + z~ (x + e,y + e 2 (y + 1)) 
+z- (x - e, y - £ 2 (y + 1)) + z (x + e, y + e 2 (y - 1)) + z + (x - e, y - e 2 (y - 1)) 
-z + (x + e, y + e 2 {y + 1)) - z + (x - e, y - e 2 {y + 1)) + z + (x + e, y + e 2 {y - 1)) ] 
+i> (x + e, y + e 2 (?/ + 1)) - if> (x - e, y - e 2 {y - 1)) + v (x - e, y - e 2 (y - 1)) 
-v (x + e, y + e 2 {y + 1)) - v (x - e, y - e 2 (y + 1)) + v (x + e, y + e 2 {y - 1)) , (1.26) 

2z+ (x, y) = 7[ - 4z + (x, y) + z+ (x - e, y - e 2 (y - 1)) + z + (x + 5, y + 5 2 (y + 1)) 
+z + (x - e, y - e 2 {y + 1)) + z + (x + e, y + e 2 {y - 1)) + z~ (x - e, y - e 2 (y - 1)) 
-z (x + e,y + e 2 (y + 1)) - z~ (x - e, y - e 2 (y + 1)) + z~ (x + e, y + e 2 (y - 1)) ] 
+2ip (x, y) - ij) (x + e, y + e 2 {y + 1)) - x/j (x - e, y - e 2 {y - 1)) 
+v (x - e, y - e 2 {y - 1)) + v (x + e, y + 5 2 (y + 1)) 

-v (x - e, y - e 2 (y + 1)) - v (x + e, y + e 2 (y - 1)) , (1.2c) 

v (x, y) = u) 2 [ - 2z + (x, y) + z + (x - e, y - e 2 {y - 1)) + z + (x + e, y + e 2 {y + 1)) 
-z- (x - e, y - e 2 (y - 1)) + z~ (x + e, y + e 2 (y + 1)) ] + y[v (x + e,y + e 2 (y+ 1)) 
+v (x - e, y - £ 2 (y - 1)) + v (x - e, y - e 2 (y + 1)) + v (x + e, y + e 2 (y - 1)) 
-4v(x,y)]. (1.2(f) 

Finally, straight forward differentiation in the continuous coordinates x and y, up to 
0(e 2 ), lead to a set of partial differential equations for the time evolution of these four 



\j = Vi+ij + Vjj_i + Vi_ij + Vij+i - 4vi,j . 



(A.l) 



+v (x + e, y + e 2 (y - 1)) - 4v (x, y) ] 



(1.2a) 
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fp = Aeu 2 z- + 2e 2 (u 2 (z+ + 2yz~) + 7 v xx ) , (1.3a) 

z = eip x + e 2 {^z~ xx + yip y ) , (1.36) 
1 
2 

v = 2ecu 2 Zx + e 2 (u 2 (z+ + 2yz,~ + 2z+) + 2 7 v ra .) . (1.3d) 



e 2 ( -^xx + lKx-^y + ^y ) > (l-3c) 



Appendix B. Derivation of the differential equations: x — 0. 

In this appendix, we derive a set of partial differential equations for the evolution of the 
diagonal covariances. The dynamics in the boundary {x = 0} is different from the bulk 
dynamics due to the stochastic collisions and is not related to the physical boundary 
conditions concerning the coupling with the heat baths. 

The dynamics along this boundary are obtained by considering the difference 
equations fl29al) - fl29a1) along the diagonal (i = j), and along the the sub-diagonal 
(i — j — 1). Along the diagonal the difference equations become 

Qi = 2^ 2 (-2s+ t + z+ M + z+_ 1 + z- +M + zr j _ 1 ) +7(r j „ 1 + r i+1 -2T J ) , (2.1a) 
s+=7(-2s+ + z+ M + z+ i _ 1 -zr +M + z - i ^ + ^^^ , (2.1b) 

f l = 2u 2 (-si + z+_ ltl + z- +1 )+ 1 (T l . 1 +T i+1 -2n . (2.1c) 

The continuum limit rule (13"3"|) for the diagonal on sub-diagonal matrix elements 
can be written as 

m i+AiA+Aj = m(/e, -1 + e 2 {-fy + s)) , (2.2) 

with the shift functions / and s defined as before. By using this rule and differentiating 
the resulting continuous equations, we obtain up to 0(e 2 ), 

fi = Aoo 2 (z + - s + + z") + Aeuj 2 (zt + z x ) + 2eV(z+ + z xx + 2y(z+ + z~)) , (2.3a) 

s + = Q - if, + 2 7 (z + - s + ) + e( - tp x + 2 7 z+) + e 2 { - -ip xx + 7 z+ 

-(y + + 2 7 (yz+ - z~) + 2v y ) , (2.36) 
t = 2^(z + -s + +z^) + 2^ 2 (z+ + z x ) + eV(z+ + z xx + 2(y + l)(z+ + z~)) . (2.3c) 

Analogously, on the lower diagonal (i — j — 1), the difference equations become 





= u 2 ( 


- 4Z+ + S+ + Z+ + S+ ! + Z+ _ 2 - 


"1" Z i+l,i-l + Z i,i-2) 


(2.4a) 




= 7( 


~ ^ Z i,i-1 + Z i+l,i-l + Z i,i~2 + S t ~ Z t+l,i 


-1 ~ S t-l,i-l + Z t,i-2) 








— &i + Ti — Vi+i^-i — Tj-i - 


f V i: i„ 2 , 


(2.46) 


2z+_i 


= 7( 


~~ 4 z i[i-i + + z^ + s^_ 1 + z i(i _ 2 - 


" Z i+l,i-l + Z i,i-2) 








+2^i,i-i - ^1+1,1-1 - + Ti + v m>i _ 


-1 — — , 


(2.4c) 


Vi,i_i = 


= - 2 ( 


" 2<i_! + s+ + z+ + z r + + 7 ( 


Vi+l,i-l + Vj 5 j_ 2 — 2Vj j j_i) , 


(2.4a 1 ) 
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and by using (12. 2p and keeping differential terms up to 0(e 2 ) we arrive at 



z/> = 2u 2 ( - z+ + s + + z") + 2b(2uj 2 z~ + 7 v x ) + e 2 (2^ 2 (z+ + 2z" + z+ - s 



+ (2y - l)z;) + ^3v xx + 2j/v„)) , (2.5a) 

z = -{if>-to) - 7z~ + ei/>x + e 2 {ip xx +y^y-w y + T y + 7(z^-z+ + s++z~)) , (2.56) 

z' + = ^(V'-^)-7(z + -s + )+^(-^--^+Vj / +T y + 7 (z++z+-s+-z;)) , (2.5c) 
v = uj 2 { - z + + s+ + z~) + 2e(u; 2 z; + 7V;c ) + £ V(z+ + 2z; z + 2z+ + 2yZy) 




The set of partial differential equations (I2.3a[) -( f273cj) and (I2.5al) -( l2~5dl) . describe the 
dynamics of the covariance on the line x = 0. 

Appendix C. Derivation of the differential equations: y — ±1 

Along the boundaries, where the system is coupled with the heat baths, the dynamics 
is different from that in the bulk not only because of the coupling with the bath itself, 
but also because the deterministic force felt by the edge oscillators feel from the bulk 
is not counter-balanced and, moreover, the stochastic collision at the edges is a "one- 
sided" process. Clearly, the dynamics depend on whether the BC are free or fixed. 

In this section we derive the covariance dynamic equations at the boundaries 
y = ±1. These boundaries correspond, on the original matrix, to the matrix elements 
of the first and last rows and columns. Restricted to the semi-infinite plane x > 0, 
i.e., to the domain ( }3"%j) . y — — 1 corresponds to the first matrix column, while y = 1 
corresponds to the last matrix row. Here we specialize on the boundary y — — 1. The 
BC at y = 1, that can be obtained analogously, lead to the same information that is 
extracted from the BC at y = — 1. 

Appendix C.l. Fixed boundary conditions: y = — 1 

We recall that fixed BC are defined by the relations Agi = q\ and A.qN+i = —qN- Our 
starting point is transforming the set of difference equations f|T5l) - f|T7j) into a set on the 
variables z ± and v. Note however, that the domain of Agj, which is i G [1, A + 1], is 
different from the domain of Pi, i G [1, A]. As a consequence, the terms Zjv+ij are non 
symmetrizable. 

In what follows, we restrict to i G [3, A — 1], where all terms are symmetrizable 
and all covariances are "far enough" from the diagonal. By finally denoting 4 1 = i/j^i, 
= zf x and u it i = v^i, we obtain 



& (1 = uj 2 [ - 4c+ + c + i,i + c+ 1,1 + cr+i,i - Ci,i + < 2 - - Ki 



+7[%i,i + ^-1,1 + v i,2 - 3i/ M ] , 



(3.1a) 



2Q = 7 [ - k~i + c+i,i + Cm + c; - c&i,i - cti,i + \2 + < 2 ] 



(3.16) 



Dynamics of anomalous heat transport 23 
2ft = 7[ " 3ft + C4i,i + C + -!,i + Ci - C m ,i " Cw + <2 + z" 2 ] 

— ^[CvL + C\l] + 2 4>i,l + 0t+l,l - ^i,2 + - + V ij2 , (3.1c) 

v itX = u 2 i- 2cj + cii,i + + zj" 2 - 17,2] - Ki 

+7^+1,1 + + v i,2 - 3v i)2 ] • (3.1(f) 

= w 2 [ - 4z+ + z+ + z+ 1>2 + z+ 12 + C+! - zr 3 + z" +li2 - zr 1)2 + Q] 

+7[v<+i,2 + Vi |3 + Vi_i i2 + i/ 4| i - 4v; j2 ] , (3.1 e) 

2Z i " 2 = 7[ - 4z l " 2 + Z *~3 + ^+1,2 + ^-1,2 + C"l + z t 3 - Z i+1,2 - Z tl,2 + Cl] 

+^i+l,2 - ^i,3 + Vi j3 - V i+ i )2 - V;_ lj2 + V ifl , (3.1/) 

2z ^2 = 7[ - 4z+ + z+ + z+ 1>2 + z+ 12 + Cj + zr 3 - z" 12 - z r 1)2 + 

+2^ ij2 - ^ i+ i >2 - ^i,3 + v i)3 + v,+i )2 - Vi_i )2 - , (3.1#) 

V i)2 =U 2 [- 2z+ 2 + Z+3 + Z+ 1)2 - zTg + Z^ 1)2 ] 

+7[ v i+i,2 + v i>3 + Vi_i )2 + 2^,1 - 4v i)2 ] • (3.1/i) 

In the continuum limit, the column index 1 corresponds to y — — 1. As a result, (1331) 
becomes 

m i+ Au+Ai = m(x + /e, -1 + £ 2 (-/ + s)) , (3.2) 
with the shift functions / and s defined as before. Using this and proceeding as in 



Appendix A , the dynamics at y = — 1 for fixed BC is described by the following set of 
partial differential equations 

= cj 2 (z + -z -2C + ) - lu - 7(1/- v) + e[- cj 2 z+ - 7 v x ] , (3.3a) 

2C" = 7 ( Z+ + Z "-C + -D - KC + C) -^ + 0-2z/ + v- £[ 7 z+-vJ , (3.35) 

2C+ = 7 (z + + z~-C + -D - ^(C + + C) - ^ + + v + e[- 7 z+ + 2i/ x -vJ , (3.3c) 

v = w 2 (z + -z--C + + C") - lv - 7(f-v) + e[u; 2 (C + - z ^)-7v,] , (3.3d) 

^ = cu 2 (-z + -z" + C + + C) + 7(^-v) + 4oj 2 z~ + 2e 2 [u; 2 z+ + 7 v x J , (3.3e) 

2z~ = 7 (-z + -z- + C + + C") + ^-v + 2^ , (3.3/) 

2z + = 7 (- z + -z- + C + + C") - + 2e 2 [ 7 z+ + 2v„] , (3.3s) 

v = 7(1/ - v) + 2^ 2 z" + £> 2 (z+ + 2z+) + 2 7V;c J , (3.3/i) 



Appendix C.2. Free boundary conditions: y = — 1 

We start by recalling that free BC are defined by setting Ag x = Ag^r +1 = 0. This means 
that in this case, the domain of the phase variables is % G [2, N] for Ag« and i G [1, N] 
for Consequently, while Zj >1;1 is well defined, its symmetric component z^j is not, 
thus restricting the symmetrization of the covariance z (1281) . At variance with the case 
of fixed BC, here it is necessary to consider non-symmetrizable terms. In analogy to the 
previous section, we distinguish boundary covariances from their bulk counterparts by 
defining = v^i and Q^i = z^i, recalling that is non-symmetrizable. Moreover, 
note that for free BC, ip has no boundary component, as the dynamics of T/>i, 2 corresponds 
to that in the bulk (see (ED). 
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By transforming the set of difference equations (j2ip - (!23p into the covariance 
variables of section [2731 and restricting to i G [4, N — 1], we obtain 

ipi,2 = uj 2 [- 4z+ 2 + z+3 + z+ 12 + z+ lj2 - z" 3 + z^ 1>2 - zr_ 12 + Ci,i] 

+7 [vi+1,2 + v, i3 + Vi-1,2 + ^,i - 4v i)2 ] , (3.4a) 

2^2 = 7 [- 4z *~2 + zr 3 + z r +li2 + z r lj2 + z + - z+ 12 - z+ 12 + 

+^i+l,2 - ^i,3 + Vi )3 - v m>2 - Vi_ 1)2 + , (3.46) 

2 <2 = 7 [- 4Z ^2 + <3 + Z i+1,2 + Z tl,2 + Z *~3 ~ Z i+1,2 ~ Z i-1,2 + Ci,l] 

+2^ 2 - - ^t,3 + v i,3 + Vi+1,2 - Vi_i i2 - v itl , (3.4c) 

Vi,2 = UJ 2 [-2Z+ + Z+ + Z+ lj2 - Z- 3 + Z m J 

+7 [vt+i.2 + v i)3 + Vi_i j2 + - 4v i>2 ] • (3.4c?) 

0,i = 7 [ z ^2 + z »\2 - Ci,i] - Ki,i ~ 4>i,2 + - Vi-i,x + v i)2 , (3.4e) 

z>i,i = cu 2 - 0,i + z ^2 - z i^] - + 7 h+i,i + ^i-i.i + v i,2 + — 3z/ i; i] . (3.4/) 

Furthermore, using (13.21) . we obtain the continuous version of the equations above. 

Differentiating these later equations and keeping terms up to 0(s 2 ), we obtain the 
following set of partial differential equations: 

ip = u 2 {- z~ + C-z + ) + 7(v-v) + 4:£u 2 z~ + 2e 2 [u 2 z+ x + -fv xx ] , (3.5a) 

2z~ = -f(-z~ + (-z + ) + v-v + 2eip x , (3.56) 

2z + = 7 ( -z- + C-z + ) - v + v + 2e 2 (u 2 zt x + 2v„), (3.5c) 

v = 7(1/- v) + 2ea; 2 z- + £V(z+. + 2z+) + 27V,.,.] , (3.5(f) 

C = - 7 (-z~ + C-z + ) -K + v-^ + e[- 7 (z+ + z") + v x - v x + fa] , (3.5e) 

v = u 2 (z + -z-) - 7(1/ - v) - lu + e[uj 2 (z- + ( x -zi) - ^ x ] . (3.5/) 

We remark that the covariance's equations at y = — 1 for the diagonal terms, namely 
those obtained for i = 1 and i = 2, yield the same solution to leading order, with the 
additional constraint = T + . 
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